The Learning Curve

home *** CD-ROM | disk | FTP | other *** search

/ The Learning Curve / The Learning Curve (Weird Science, 1996).iso / electronics / analog_emulator / aces-docs / appendixd < prev next >

Wrap

Text File | 1995-09-07 | 40.3 KB | 1,336 lines

** APPENDIX D ** (Preliminary Docs.) -Example 01- -SpringMass- This classical example can best illustrate analog programming techniques. Consider the second-order differential equation describing tne motion of a mass (M) which is connected to a stationary support by the combination of a spring of stiffness (k) and a damping device with viscous damping factor (b). Using the basic equation F = MA and writing a force balance we have: (b Xdot + k X) - (M * Xdotdot) = 0. where: Xdotdot = acceleration of mass. (ft./sec./sec.) Xdot = velocity of mass. (ft./sec.) X = deflection of mass from rest position. (ft.) M = units of mass. (lb.- weight) k = spring constant. (lb./ft.) b = damping factor. (lb.-sec/ft.) Rewriting in terms of the highest derivative since we will use integration.. Xdotdot = (Xdot * s/-M) + (X * k/-M) or Xdotdot = -1/M (s Xdot + k X) Assuming a mass of size 5.0 (161.0 lb. weight) For initial conditions, we have: (X t=0) = 0.0 ft. / sec. (at rest). (X t=0) = -10.0 ft. (initial deflection). For the constants let's use: k = 1 lb. / ft. b = 2 lbs./ sec./ sec. Assumed maximums: X = 10.0 ft. / sec. X = 10.0 ft. The analog scaling becomes: [Xdotdot/10] = -1/M (b [Xdot/10] + k [X/10]) With scaled initial conditions of: ([Xdot/10] = 0.0, t = 0 ) {at rest} ([X/10] = -1.0, t = 0) {-10 / 10} Independent value (time units) will be in seconds. Since the undamped frequency will be k/m (radians/second) and the total length of time to observe 2 cycles would be approximately 25 sec.. For the step size, we can use .05 sec.(total time/ 20 steps for each cycle). We can accept a print frequency of 2 seconds. Preparing the block diagram: We start by implementing the equation for the highest derivative. We assume we have the input terms of Xdot and X available. We can use a weighted summer (W) block to form the X double dot variable. We can assign Block 1 to this module. Input 1 will come from a block representing X dot and Input 2 will come from a block representing X. Parameter 1 is associated with Input 1 and Parameter 2 is associated with Input 2. (Parameter 3 and Input 3 are not used.) Drawing the diagram at this point we have: p1 e1 --O-|\ 1 | \______ [Xdotdot/10] | / e2 --O-|/ [W] p2 Notice that the output of Block 1 is labeled with the scaled quantity X double dot. Having formed the output of X double dot we can now use it to form X dot. X dot is formed by integrating X double dot. We can use an Integrator module for this operation. Since the parameter used is 1.0, we can use the first input of the Integrator module (I). We can assign this module Block 2. We draw a line from the output of Block 1 to the first input of Block 2. The output of Block 2 will be labeled with the scaled quantity of X dot. The first parameter (Parameter 1) is used to set the initial value of X dot. p1/IC = 0.0 \ O 1 >----|-|\ / |-| \ 2 p2 O|-| \______ [Xdot/10] [Xdotdot/10] |-| / |-| / [I] p3 O|-|/ We now can form the quantity X by integrating the output of X dot. We do this with another Integrator module, which we will assign to Block 3. Again we have a parameter of 1.0 associated with the scaling so we can use the first input of Block 3. The output of Block 3 now represents the deflection X. The first parameter (Parameter 1) is used to set the initial value of X. p1/IC = -1.0 \ O 2 >----|-|\ / |-| \ 3 p2 O|-| \______ [X/10] [Xdot/10] |-| / |-| / [I] p3 O|-|/ We now go back where we started. We assumed we had the quantities X dot and X when we formed X double dot. we now have the actual block assignments for these terms. We now connect a line between the output of Block 2 (X dot) and Input 1 of Block 1. We also connect a line between the output of Block 3 (X) and Input 2 of Block 1. Parameters 1 and 2 will contain the respective coefficients for -b/m and k/m terms. The scaled value of these parameters can be written on the diagram for reference. p1 = -(b/M) [Xdot/10] \ 2 >--O-|\ 1 / | \______ [Xdotdot/10] \ | / [X/10] 3 >------O-|/ [W] / p2 = -(k/M) Our configuration diagram is now complete. ----------- |\ /| | \ / | | \ / | | \ / | | X | Figure X | / \ | | / \ | | / \ | |/ \| ----------- We now have all the information, at hand, that we need to prepare an ACES dataset to implement the above scaled equations and run the analog. The first (configuration) data section will become: acc.[Xdot2/10] 1 W 2 3 0 veloc.[Xdot/1] 2 I 1 0 0 distance [X/1] I 3 2 0 0 # The second (parameter) data section will become: b/m,k/m 1 1.0 2.0 0.0 10.0 0.0 Xdot(t=0) 2 0.0 0.0 0.0 10.0 0.0 X(t=0) 3 -1.0 0.0 0.0 10.0 0.0 # Integration data section: 0.05 25.0 0.2 # Print assignment data: 1 2 3 0 0 0 0 # Plot assignment data: 76 1 2 3 0 0 # Final record mark (include 1 extra mark for end function data area): # This is the minimum data required to run an ACES application. The plot data can be expanded and arbitrary function data added if required. (The dataset above is named SpringMass and is included in the Examples drawer of your distribution). -EXAMPLE 02- -SineWave- This example is used to show the effects of chosen step size on simulated analog accuracy. ACES includes both sine (S) and cosine (C) modules. Here the sine-omega-t and cosine-omega-t functions are generated analyticly and then compared against the output of of the included functions. All ACES modules use the AMIGA double-precision floating math library functions. The equations for this example are: Sine-omega-t = (integral) Cosine-omega-t Cosine-omega-t = -(integral) Sine-omega-t Since the values of sine and cosine are in the range -1.0 - +1.0 The scale factor will be 1.0... Implementing the analog flow diagram we have: ----------- |\ /| | \ / | | \ / | | \ / | | X | Figure X | / \ | | / \ | | / \ | |/ \| ----------- Adding the internal math functions and subtracting to form error term: ----------- |\ /| | \ / | | \ / | | \ / | | X | Figure X | / \ | | / \ | | / \ | |/ \| ----------- We will check at 1 hertz (1 cycle/second = 6.2832 radians/sec.). Preparing a dataset: A dataset named SineWave has already been prepared. SineWave can be found in the Examples drawer. The end value of t = 3 sec.. is chosen so that three cycles will be calculated. The step size chosen (0.05) will give 20 steps/cycle or 60 steps for the run. This is a reasonable step size and should give good results in terms of both frequency and amplitude accuracy. Variables will print out every 0.1 second of operate time. All 6 variables are printed. The plot file has been set up for variable presentation and scaling. Load and run the SineWave dataset. Select Operate Mode and observe the calculated data and especially the error calculations. Return to the Setup mode and select EDIT from the DATA menu. Click on the Time Parameters box and change the line to change the step size to 0.1 (10 steps/cycle) as : 0.1 3.0 0.1 Click on Modify to change the data. Exit the editor and re-Run the dataset. Select the Operate mode and compare results to the previous calculation. Return to Setup mode, select the editor, and change the step size to 0.025 (40 steps/cycle) as: 0.025 3.0 0.1 As before, Modify, Exit, ReRun, and Re-Operate. Compare the results of the the calculations to the two previous times. Get a feel for the physical time involved for doing the calculations. Compute time and accuracy are somewhat inversely proportional to each other and are a function of step size. At an artificially large step size, the calculations will fall apart primarily due to the phase errors. Return to the Setup mode and select PLOTTER. Set up for two channels of Y, and select the highest plot Density. Put the Plot window to the rear of the window and reRun the dataset. Select Operate mode and then put the plot window to the front. You should see three cycles of the sine and cosine function being plotted in the Plot window. The computation will proceed to the end and reach the Hold mode. Put the plot window to the rear. Select Setup mode, Enable the Editor and change the first and second lines in the plot record to: 2 1 2 0 0 0 1.0 -1.0 1.0 -1.0 Close the Editor, reRun the data, and select Operate mode. Bring the plot window to the front and observe the plot. We are now plotting the sine function verses the cosine function (a circle). Get familiar with the plotter, change scaling, add grids, change density etc. -Example 03- -PHYSBE- The differential equations in this example were first introduced by McLoed as a PHYsiological Simulation Benchmark Experiment (PHYSBE) for comparing simulation tools and techniques available for physiological calculations. The equations represent a very simplified volumetric model of the human circulatory system including the heart's pumping function. This example illustrates the use of ramps and arbitrary functions in analog computation. PHYSBE has been used as a de-facto benchmark to test analog equipment and particularly analog emulators for speed and accuracy. This might be considered a "medium" scale simulation, based on the number of modules (56) and complexity of the system. A continuous material balance is performed during computation. The differential equations to be solved are as follows: Heart right ventrical- dVrv ----- = Ftv - Fpv Vrv (t=0) = 150.0 cc. dt Prv = Vrv/Crv Pulmonary arterial system- dVap ---- = Fpv - Fps Vap (t=0) = 120.0 cc. dt Pap = 0.133 Vap Fpv = 90.0 (Pav - Pap) Pulmonary venous system- dVvp ---- = Fps - Fmv Vvp (t=0) = 240.0 cc. dt Pvp = 0.033 Vvp Fps = 0.70 (Pap - Pvp) Heart left ventrical- dVev ---- = Flv -Fmv Vlv (t=0) = 150.0 cc. dt Plv = Vlv/Clv where Clv is a function of time. Flv = 17.0 (Pvp -Plv) Aorta- dVao ---- = Fav - Fas Vao (t=0) = 150.0 cc. dt Pav = 0.80 Vao Fav = 80.0 (Plv -Pav) Systemic circulation- dVsc ---- = Fas - Fvs Vsc (t=0) = 3340.0 cc. dt Psc = 0.0153 Vsc Fas = 1.630 (pav - Psc) Vena cava- dVvc ---- = Fxx - Fxx Vvc (t=0) = 500.0 cc. dt Pxx = 0.0040 Vvc Fvs = 1.650 (Psc - Pvc) Where all volumes are in cubic centimeters, all pressures in mm. hg., and all flows are in cubic centimeters / second. We will add a material balance to sum all the blood volumes considered. Vt - (Vrv + Vap + Vvp + Vlv + Vao + Vsc + Vvc) = 0.0 The scaled equations become: Heart right ventrical- [dVrv ------- = [Ftv/250] - [Fpv/250] [Vrv/250] (t=0) = 150./250 = 0.6 dt/250] [Prv/250] = [Vrv/250]/[Crv/250] Pulmonary arterial system- [dVap ------- = [Fpv/250] - Fps[ Vap (t=0) = 120./250. = 0.48 dt/250] [Pap/250] = 0.133 [Vap/250] (Pap limited to positive) [Fpv/250] = 90.0 ([Pav/250] - [Pap/250]) Pulmonary venous system- [dVvp ------- = 0.5 (Fps/250] - [Fmn/250]) Vvp (t=0) = 240./500. = 0.48 dt/500] [Pvp/250] = 0.033 ( 2.0 [Vvp/500]) = 0.066 [Fps/250] = 0.70 ([Pap/250] - [Pvp/250]) Heart left ventrical- [dVev ------- = [Flv/250] - [Fmv/250] [Vlv/250] (t=0) = 150./250. dt/250] [Plv/250] = [Vlv/250]/[Clv/250] where Clv is a function of time. [Flv/250] = 17.0 ([Pvp/250] -[Plv/250]) Flv limited to positive. Aorta- [dVao ------- = [Fav/250] - [Fas/250] [Vao/250] (t=0) = 150./250. = 0.6 dt/250] [Pav/250] = 0.80 [Vao/250] [Fav/250] = 80.0 ([Plv/250] - [Pav/250]) Systemic circulation- [dVsc -------- = 0.05 ([Fas/250] - [Fvs/250] dt/5000] [Vsc/5000] (t=0) = 3340./5000. = 0.668 [Psc/ = 0.0153 (5000/250) [Vsc/250] = 0.306 [Vsc/250] [Fas/250] = 1.630 ([Pav/250] - Psc[250]) Vena cava- [dVvc -------- = 0.25 [Fvs/250]- [Ftv/250] [Vvc/1000] (t=0) = 500./1000. = dt/1000] 0.5 [Pvc/250] = 0.0040 (1000/250) [Vvc/1000] = 0.016 [Vvc/1000] [Fvs/250] = 1.650 ([Psc/250] - [Pvc/250]) Material balance- [Vt/5000] - ( 0.05 [Vrv/250] + 0.05 [Vap/250] + 0.10 [Vvp/500] + 0.05 [Vlv/250] + 0.05 [Vao/250] + 1.00 [Vsc/5000] + 0.2 [Vvc/1000]) = 0.0 The scaled analog diagram is: ----------- |\ /| | \ / | | \ / | | \ / | | X | Figure X | / \ | | / \ | | / \ | |/ \| ----------- The Heart function Cv's have been implemented on two function modules driven by a ramp function. Type 1 curve contains the right ventricle function and Type 2 contains the left ventricle function. The implementation is shown in the following figure: o---------o |\ /| | \ / | | \ / | | \ / | | X | Figure X | / \ | | / \ | | / \ | |/ \| o---------o The total run time will be equivalent to 4 beats of the heart or 4.0 seconds. A step size of 0.01 sec has been deemed appropriate and a print interval of 0.1 seconds will suffice. A dataset named PHYSBE is in the Examples drawer. Note that the solution follows a stable cycle condition and does not "drift" due to small errors and unbalances as it probably would on an actual analog machine. -Example 04- -AnalogPID- ACES has a Proportional-Integral-Derivative (PID) analog controller-type module in it's complement of module types. The PID module (%) is implemented, internally, as a macro of ACES standard modules. The following example uses one of the three available % modules in a simple control scheme that maintains the height of liquid in a tank. A step change is made in the flow rate out of the tank and the controller must change the input flow rate to maintain level. The level control is complicated by the fact that there is a measurement delay associated with the tank liquid level. The controller is to be adjusted to minimize overshoot and settling time due to the outflow step change. This example demonstrates use of the PID (%) and the Unit Delay (U) modules. The tank liquid level is to be maintained at 5.0 feet. The outlet flow rate will be increased by a step change of 25% one second after the start of the computation. Only Proportional and Reset controller modes will be utilized. The equations for this example are: Qout (t=0) = 0.1 H/sec. Qin = 0.2 Zout H/sec. Qout (t=>1) = 0.125 H/sec. dH --- = (Qin - Qout) H (t=0) = 5.0 feet. dt H = f(Qin-Qout,t) Zin = H/10.0 = 0 < 1.0 Zout = 0 < 1.0 * Q is in units of tank height. (H/sec.) Z refers to the controller which is assumed to have a maximum input/output of 1.0 . The scaled equations are: [Qout/10] (t=0) = 1./10 = 0.1 [Qin/0.2] = Zout [Qout/10] (t=>1) = 1.25/10 = 0.125 [dH ------ = ([Qin/10] - [Qout/10]) [H/10] (t=0) = 5./10. = 0.5 dt/10] Initial Controller settings: Kp = -10.0 (Inverse action) Ki = 2.0 repeats/sec. Kd = 0.0 (No action) Zout (t=0) = 0.5 Zsp = 0.5 (Equivalent to 5.0 ft. level) A total operation time of 10 seconds is to be observed, with the outflow increased by 25%, 1 second after starting operation. The measurement delay is set at 1.0 second. The analog flow diagram is as follows: ----------- |\ /| | \ / | | \ / | | \ / | | X | Figure X | / \ | | / \ | | / \ | |/ \| ----------- A dataset named AnalogPid can be found in the Examples drawer. After "Running" the data set, select the Control Panel option. The PID controller settings can be varied for each run by changing Par1 and Par2 (proportional and integral gains respectively) and re-operating to get a new condition. Rep/Op mode may also be used to cycle the example automatically between I/C and Operate modes. Select Hold mode to halt Rep/Op. Experiment with the controller parameters and try to optimize the settings to reduce tank liquid level overshoot and settling times. Note: Parameter settings made with the Control Panel are only temporary. To make any parameter changes permanent, you must Edit and Save the dataset. -Example 05- -SamplePid- This example uses the same simple liquid tank level problem as in the previous example 04. This example uses sample-hold (zero-order-hold) and unit delay hybrid analog devices to implement the controller. The implementation of the sampling controller is as follows: o---------o |\ /| | \ / | | \ / | | \ / | | X | Figure X | / \ | | / \ | | / \ | |/ \| o---------o The control algorithm is implemented by using zero-order-hold (Z) and unit delay (U) modules. Parameters Ka and Kb are somewhat equivalent to their proportional and integral gain analog PID counterparts. The analog flow diagram is as follows: o---------o |\ /| | \ / | | \ / | | \ / | | X | Figure X | / \ | | / \ | | / \ | |/ \| o---------o The dataset for this example can be found in the Examples drawer under the name SamplePid. For all practical purposes, this example can be operated very similarly to Example 04. -Example 06- -Splat- The pilot ejection problem, featured in this example, has been used as a basis for comparison of analog simulations using continuous system languages such as MIMIC, CSSL, MIDAS, and others. The problem statement is taken from the MIMIC programming language manual. In this investigation, the trajectory of a pilot and his seat assembly, while ejecting from the aircraft, is calculated. It is necessary to determine whether the pilot is safely ejected clear of the trailing part of the aircraft without subjecting him to unnecessary forces. The vertical stabilizer is located 14 feet higher and 60 feet to the rear of the cockpit. The seating system is mounted on rails and is designed to move the pilot and his seat assembly, at an angle, up and towards the rear of the craft, at a specified velocity. The pilot and seat assembly clear the aircraft after traveling a vertical distance of 4.0 feet along the rails. At this point, the pilot is subjected to the drag forces due to his velocity in the air and follows a trajectory out and, hopefully, above the path of the plane. The equations involved in this example are: dTheta ------ = 0.0 (y < 4) Theta (t=0) = -.04532 rad. dt (from vertical) = -32.2(cos Theta)/V (y > 4) dX -- = V (cos Theta) -900. X (t=0) = 0.0 fps. dt dY -- = V (sin Theta) Y (t=0) = 0.0 fps. dy dV -- = 0.0 (Y < 4) X (t=0) = 900.0 - Vs dt = -0.001698 * Vsqr -32.2 (sin Theta) (Y > 4) We will include Yp, a function of X of the planes upper surface, to aid in visualization of the pilot's trajectory on the plot. Yp = f(X) The scaled equations are: [dTheta ------ = 0.0 (y < 4) Theta (t=0) = .04532 rad. dt/1] = -32.2/1000 (cos Theta)/[V/1000] (y > 4) [dX ------- = (1000/100) [V/1000] (cos Theta) -900/100. dt/100] X (t=0) = 0.0 [dY ------ = (1000/20) [V/1000] (sin Theta) dy/20] [Y/20](t=0) = 0.0 [dV -------- = 0.0 (Y < 4) [V/1000](t=0) = 900./1000 - [Vs/1000] dt/1000] = -0.001698 * (1.0e03)([Vsqr/1.0e06] (-32.2/1.0e6) (sin Theta) ([Y/20] > 4) [Yp/20] = f([X/100]) Escape velocity and rail inclination are the parameters to be manipulated. The object of this study is to find the minimum velocity and the best angle to use. The aircraft is assumed to be in level flight at an airspeed of 900.0 FPS. Use the Control Panel to manipulate the parameters and the QUIK-PLOT window to observe trajectory. The ordinate represents -X (the distance along the plane) and the abscissa Y (height above the pilots normal seating position). Use of the logical AND (&) and the quit module (Q) module are demonstrated. The analog flow diagram is as follows: o---------o |\ /| | \ / | | \ / | | \ / | | X | Figure X | / \ | | / \ | | / \ | |/ \| o---------o A dataset named Splat is in the Examples drawer. -EXAMPLE 7- -Plate- In order to solve partial differential equations on an analog computer, it is generally found that the best method of approaching the solution is to reduce the partial equations to a set of ordinary differential equations. The one-dimensional diffusion equation commonly found in heat transfer work can be broken down to one continuous variable and one space variable with time being the independent variable. We have the following partial differential equation describing the one-dimensional diffusion of heat in a solid: (partial)sqd. T k (partial) T --------------- = ------------- (partial) Xsqd. (partial) t Dividing X (the thickness of the plate) up into n equal spaces yields the following equations for heat flow (through the thickness of the plate) at the n'th station: dT(n) (Qin - Qout) * A ----- = ----------------- dt rho * cp * V(n) k * A Qin - Qout = -------- * ( T(n-1) - 2T(n) + T(n+1) delta X Where k = conductivity in BTU / lb. - ft. - deg. F. rho = density in lb. / cu. ft. cp = heat capacity in BTU / lb. - deg. F. A = area considered in sq. ft. V = delta X * area of each section. delta X = internode distance (section size) Q = BTU / hr. sigma = Stefan-Boltzman constant 0.174e-08 (radiation) Equations for the internal sections: dT(n) (Qin - Qout) * A ----- = ------------------ dt rho * cp * V(n)/2 k * A Qin - Qout = ----------- * ( T(n-1) - 2T(n) + T(n+1) delta X Equations for the surface section: dT(0) (Qin - Qout) * A ----- = ----------------- dt rho * cp * 1/2 V(n) k * A Qin - Qout = ----------- * ( T(1) - 2T(0) + T(1) + 1/2 delta X k * [Ta - T0 + epsilon * E * ((Ta + 460) - (T0 + 460)) For the symmetrical center station: dT(5) (Qin - Qout) * A ----- = ------------------ dt rho * cp * V(n)/2 k * A Qin - Qout = ----------- * ( T(4) - 2T(5) + T(4) delta X / 2 Consider the following assumptions: Heat conduction is equal across the area of the plate (one dimensional). We can use symmetry and calculate either half of the plate (Both surfaces are subject to the same conditions of heat flow by convection and radiation). The material consists of a mild steel plate 1.0 inches thick radiating to a semi-infinate heat absorbing structure (constant temperature). The thermal coefficients are constant over the temperature range. The temperature gradients are "mild" and can be approximated by straight lines, without excessive error, across the number of sections considered. We will use 5 sections to represent the plate half width with 5 temperature nodes. The internode distance would be 1 inch / 2.0 divided by 4 or 0.125 in. across each section. (The surface and center sections are half sections/volumes). It is assumed that the AISI 1010 carbon steel plate, which was heated to a uniform temperature of 1500.0 degrees F., will be removed from a furnace and placed on edge for cooling. The following heat transfer coefficients are to be used: rho (density) = 490.0 lb./ cu. ft. cp (heat capacity) = 0.16 BTU / lb. - deg. F. k (conductivity) = 16.0 BTU / hr. - ft. deg. F. h (convection air) = 2.0 BTU / hr. - deg. F. E (emissivity ) = 0.80 We will scale an internal section: [ dT(n) 2.0e04 [(Qin - Qout)/2.0e04] * 1 -------- = ------- * -------------------------------- dt/2000] 2000. 490. * 0.12 * 0.125 * 0.125 /144 [Qin -Qout 2000 12.0 * 1 [T(n-1) -2T(n) + T(n+1) ---------- = ------------- * -------- * ----------------------- 2.0e04] 2.0e04 0.125/12 2000 ] Scaling the surface section: [ dT(0) 2.0e04 [(Qin - Qout)/2.0e04] * 1 -------- = ------- * -------------------------------- dt/2000] 2000. 490. * 0.12 * 0.125 * 0.125 /144 [Qin -Qout 2000 12.0 * 1 [T(1) -2T(0) + T(1) ---------- = ------------- * -------- * ------------------- + 2.0e05] 2.0e05 0.0625/12 2000 ] 2000 [Ta - T0 ------ * 2.0 * --------- + 2.0e05 2000] 4 4 4 (2500) [(Ta + 460) - (T0 + 460) -------- * 0.8 * .174e-08 * ----------------------- 2.0e05 4 2500 ] For the symmetrical half (center) section: dT(5) (Qin - Qout) * A ----- = ------------------ dt rho * cp * V(n)/2 k * A Qin - Qout = ----------- * ( T(4) - 2T(5) + T(4) ) delta X / 2 The scaled analog flow diagram is as follows: o---------o |\ /| | \ / | | \ / | | \ / | | X | Figure X | / \ | | / \ | | / \ | |/ \| o---------o A dataset named Plate can be found in the Examples drawer. Select the QUIK-PLOT window and enable all 5 Y axis channels. Channel 1 represents the surface temperature node, 2 - 4 the internal nodes, in the physical center of each section through the plate, and Channel 5 represents the node at the center of the plate. The temperature points will be plotted against time, in minutes, on the abscicca. Run the dataset and select Operate mode. Notice how the internal temperatures lag behind the surface temperature as cooling occurs. Try changing the initial temperatures and/or some of the heat flow parameters. Remember that the temperature gradient for the model is assumed to be a straight line through each section. The model must be divided into smaller sections to properly represent thermal gradients due to more severe conditions of heating or cooling. Other surface boundry conditions can be added to this basic model. -Example 08- -Algebra- Occasionally, modeling of a system of equations, that contains algebraic equations, causes the computer to be forced to calculate an "algebraic loop". This is because a particular variable appears on both sides of an equation and must be solved for in an explicit manner. Sometimes it is possible, by manipulating the equation (by substitution or other method), to avoid solving the equations explicitly. This is the best way to handle the situation where this is possible. This example assumes this is not possible and attempts to solve the equation directly. The example chosen does, in this case, lend itself to manipulation, thus providing a check solution for the explicit case. Persons, that are familiar with actual analog computing, know that you can generally program some algebraic loops directly by introducing a small, 1st order, lag in one or more of the amplifiers in the loop. This is accomplished by adding a small capacitance in the amplifier feedback path. Noise is always a by-product of calculating algebraic loops, in this manner, and must be dealt with or tolerated when measuring outputs of these variables. ACES can't handle loops, in this manner, but we can force a similar delay action through use of a combination of the "wye" (Y) and the "vacuous" (V) elements. Converging of the loop is dependent on sign, gain, and (on digital emulation methods) step size. Positive algebraic loops, with gains greater than one, are always unstable. For this example we have a simple system of equations involving X and Y: dX dY --- = ---- + XY - X dt dt and dY dX ---- = - ---- + Y dt dt Here dX/dt depends on dY/dt and vice-versa. Let: Xmax = 10. Ymax = 10. Scaling the Equations we have: [dX [dY (10) [X [Y [X ------ = ----- + ---- --- --- - --- dt/10] dt/10] 1 10] 10] 10] and [dY [dX ------ = - ------ + [Y/10] dt/10] dt/10] o---------o |\ /| | \ / | | \ / | | \ / | | X | Figure X | / \ | | / \ | | / \ | |/ \| o---------o When diagrammed the configuration appears as above. Note the heavy lines tracing the algebraic loop. This loop will be computationally unstable in its present form. We can modify this diagram, to break the loop, by inserting the Y and V elements, in combination, to simulate a slight computational delay. ( That is, a simulated, unity gain, amplifier with a small first-order delay). This will change the configuration diagram to: o---------o |\ /| | \ / | | \ / | | \ / | | X | Figure X | / \ | | / \ | | / \ | |/ \| o---------o For comparison purposes we are going to change the equation for dX/dt by substitution: dX (dX --- = - --- + Y) + XY - X dt dt or dX --- = 0.5 (XY -X + Y) dt Scaling these equations: [dX ((10 [X [Y [X [Y ------= 0.5 --- --- --- - --- + --- dt/10] 1) 10] 10] 10] 10]) [dY [dX ------ = - ------ + [Y/10] (as before). dt/10] dt/10] The analog diagram for this system of equations becomes: o---------o |\ /| | \ / | | \ / | | \ / | | X | Figure X | / \ | | / \ | | / \ | |/ \| o---------o We will solve both sets of equations simultaneously and form an error term for comparison. The algebraic loop calculation will have a small error, the magnitude of which, is adjustable by changing the Y element error coefficient (P2). Use the dataset named Algebra in the Examples drawer. Load and run the dataset. Select Operate mode and observe, in the log, the operation of both implementations. The settings for the example are somewhat optimized and the error in the loop calculation minimal. Try increasing the allowable error in the Y element (P2) and note the action on the algebraic loop calculation. Vary the convergence factor (P3) on the Y element and note any change in the physical calculating time. P3 can be varied from about 0.1 to 0.9. Lower values gives more damping. Higher values more overshoot during convergence to error value. Values around 0.5 will generally give acceptable results. The physical calculating time can be minimized by judicious setting of P3. Vary the initial values of X and Y and the range of operate total time (t). -Example 09- Curves 1,2, & 3 Sometimes, in the course of an investigation, you may find a set of simultaneous differential equations in which the parameters are not known. Generally the values of these parameters can only be found by trial and error. One possible way is to try to vary all unknown parameters incrementally, holding the others constant, until the proper solution is found. Parameter "sweeps" of this nature are very time consuming. ACES has two things in it's repertoire that can aid in finding unknown parameters and/or boundry values. The first is REP-OP mode. Placing the emulator in repeated operation, allows the computer to continuously repeat the calculations while the operator can "twittle" the parameters of interest. The second is the iterative element (?). This element can execute a parameter sweep or act as an iterative control element to adjust a variable and/or parameter based on an algorithm producing an error term. The iterative (?) element is only functional in REP-OP mode. The example used, to illustrate the use of these features, is a basic curve fitting problem. Some data has been collected in a field study. The researcher believes that the data function Y = f(X) can be represented by a mathematic model involving a very simple differential equation. The equation is as follows: dY -- = -aY + b dX The parameters (a) and (b) are to be adjusted to best fit the collected data for 0 < Y < 10000. Scaling the equation gives: [dY [ Y b --------- = -a ------ + ----- dX/10000] 10000] 10000 The first example (Curve1) is a manual implementation to be used in REP-OP mode, in conjunction with, the Control Panel to adjust observable best fit of the equation to the data points. The data points are programmed as an arbitrary function curve. An error term is formed by subtracting the collected experimental data values from the Y term of the analytical solution. The acceptability of fit is left up to the operator's judgment. There is a performance indicator programmed to give a running R.M.S. error value curve for the entire range considered. Parameter (-a) is Px on Block x, Parameter (b) is Px on Block x. The following configuration diagram represents the above equations: o---------o |\ /| | \ / | | \ / | | \ / | | X | Figure X | / \ | | / \ | | / \ | |/ \| o---------o The second example (Curve2) automates the parameter search for parameter (a). The ? element replaces the K element representing parameter (a) from the previous example. Load, Run and Repetitively Operate dataset Curve2. Notice the interaction of parameter (a) when changing parameter (b). Notice that the overall R.M.S. error is most proportional to parameter (b) while the maximum error along the length of the curve is more sensitive to parameter (a). The following configuration diagram represents the changes: o---------o |\ /| | \ / | | \ / | | \ / | | X | Figure X | / \ | | / \ | | / \ | |/ \| o---------o The third example (Curve3) further automates the search for the best fit to the data. A pair of zero-order-hold (Z) elements are programmed to store the error values near the first (X = 0.1) and last (X = 0.6) portions of the error curve. An iterative element is added for adjusting parameter (b). Action of the ? elements is such that parameters (a) and (b) will simultaneously converge toward a point of minimum error at the value of X selected. The ? elements are programmed for a "ballpark" initial condition and a nominal response gain. The operator can stop the calculations (Selecting Hold mode) when no further improvement is observed. It is found that this model does fit the data extremely well once the proper parameters are found. Try adjusting the response gain of the ? elements. The points of error comparison (values of X) can also be changed from 0.1 and 0.6, possibly to find a point of better sensitivity to instantaneous and overall fit errors. The configuration diagram for this variation is as follows: o---------o |\ /| | \ / | | \ / | | \ / | | X | Figure X | / \ | | / \ | | / \ | |/ \| o---------o ------------------------------------------------------------------------- [end]